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Abstract 

The low temperature relaxation of the magnetization in molecular magnetic solids such as Feg is 
studied using Monte Carlo simulations. A set of rate equations is then developed to understand the 
simulations, and the results are compared. The simulations show that the magnetization of an in- 
tially saturated sample deviates as a square-root in time at short times, as observed experimentally, 
and this law is derived from the rate equations analytically. 
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I. INTRODUCTION 



The low temperature relaxation of the magnetization of magnetic molecular solids such 
as Feg has proven difficult to understand ever since the earliest experimental studies [l-3|. 
The time dependence of this relaxation is highly non-exponential, and fits to forms such 
as stretched exponentials have provided no insight even when the fits seem to be good. A 
second puzzling feature is that for short times, the relaxation is observed to follow a square- 
root behavior with time in a large number of protocols: demagnetization, magnetization, 
and hole-digging in which the magnetic field is abruptly changed after the magnetization has 
been allowed to come to an equilibrium or quasi-equilibrium state in response to a previous 
value of the applied^agnetic field. A good review of the subject is given by Gatteschi, 
Sessoli, and Villain [4,]. These authors give many more references to experimental studies 
[l 3,[5-8|, theoretical analyses [9], and Monte Carlo simulations 
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The fundamental microscopic mechanism by which the spin of an individual molecule 
changes at low temperatures (say below 50 mK) is incoherent tunneling between the lowest 
energy states. In both Feg and Mni2, the anisotropy of the molecule is of the Ising type, and 
the lowest energy states have Zeeman quantum numbers m = ±S, where S is the spin of the 
molecule. The tunnel splitting between these states is of order 100 Hz (in frequency units) 
for Feg and unobservably small for Mni2. It must be stressed that in the solid, the tunneling 
is not of the coherent fiip-fiop type seen in the NH3 molecule, and previous authors have 
examined various decoherence processes by which the tunneling dynamics of a single molecule 



change from coherent to incoherent 12|, ll3| . This is not enough to explain the observed non- 



exponential time behavior, for if there were a single characteristic time scale for relaxation 
of a single molecule, and all molecules relaxed independently, the magnetization would relax 
essentially exponentially in time with the same time scale as for one molecule. Thus, the 
non-exponential time behavior is a strong indicator that the molecules in the solid do not 
relax independently of each other. The biggest and most obvious coupling between molecules 
is the dipole-dipole interaction, and while this has been considered by many previous authors 
0, 3, Is] a complete theory is still lacking. In particular, the form has been previously 
explained by Prokofeev and Stamp Q], but as noted by Gatteschi, Sessoli, and Villain {4] it 
is unclear if it applies to all situations. These latter authors also give a heuristic argument 
for the y/i law for the particular case of the demagnetization problem. We comment further 
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or this below. 

In this paper, we report on our attempt to solve this problem. Our first approach is Monte 



Carlo simulation. In this we follow in the footsteps of Refs. [10|, lUj, and many aspects of 
our simulation and the results are very similar to those found by these authors. We then 
try to understand our Monte Carlo results by developing a set of rate equations. These 
rate equations entail the distribution of dipole fields at the molecular sites. For the specific 
problem of demagnetization, we can construct an approximate model for this distribution, 
which then enables us to solve the rate equations numerically. We find that the solution 
to the rate equations matches the Monte Carlo results quite closely. Furthermore, we can 
show analytically that that the solution obeys a square-root behavior with time at short 
times. We emphasize that as in Ref. J] we have only studied the demagnetization problem. 
Further, the scaling behavior that we find for ancillary quantities also agrees entirely with 
Ref. ^|. Thus, we can claim no priority for this result. 

We also emphasize that our model for the dipole field distributions and the rate equations 
requires no further ingredients or fitting parameters beyond those involved in specifying 
the Monte Carlo process. In this paper we only look at the problem of demagnetization 
of a spherical sample with a cubic lattice in order to minimize the complications from 
demagnetizing fields, and focus on the shape independent aspects of the problem, but we 
believe that our rate equation approach offers a method to attack a much wider class of 
problems, and in the future we hope to study other experimental protocols, sample shapes, 
and lattice types. 

The plan of the paper is as follows. In Sec. [Ill we describe the basic physical model 
underlying the relaxation jo], [3]. We then describe our Monte Carlo simulations and results 
in Sec. Illli The theory for the rate equations and the bias distribution are developed in 
Sees. [IV]and|Vl Finally, in Sec. |VT]we present our analytical solution to the rate equations, 
and the -y/t law. 



II. PHYSICAL MODEL FOR RELAXATION 



As shown in Ref. ^3[, the fundamental process that governs the dynamical behaviour of 
the spins is as follows. In a short time interval dt, the spin of the ith molecule flips from 
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m = —5* to m = S, or m = S to m = —S, with a probability 

Pflip.i = ^idt, (2.1) 

with 

Here, A is defined via the statement that i A/2 is the quantum mechanical amplitude per unit 
time for a spin to tunnel between the m = ±S states, W ~ lOEdn, where Edn is the energy 
of dipole-dipole interaction between the molecular electronic spin and the nuclear spins of 
nearby nonmagnetic atoms such as N and H which are always present in the molecules 
studied, and Ei is the energy of the m = S state relative to the m = —S state due to the 

n 

net magnetic field seen by the ith molecule. We shall refer to Ei as the bias on site i |14j |. 
For Fcg, A ~ 10^^ K, E^n ~ 1 niK, and Ei ~ 0.1 K in temperature units. (We shall set h 
and ks to unity in all working formulas, so temperature, energy, and frequency all have the 
same units.) 

The dominant feature in Eq. (12.21) is the exponential suppression of the fiip rate with the 
square of the bias energy Ei, and a large part of this energy arises from the dipolar field 
of the other molecular spins in the solid, which can be estimated to be of order 100 Oe for 
near neighbour spins, leading to the energy scale 0.1 K quoted above. More explicitly, the 
dipolar part of Ei is given by 

Ei,dip = ^Kijaj, (2.3) 

(2.4) 

' ij \ ij / 

Here, Edm is the energy scale of interaction for near neighbours, a is the near-neighbour 
distance, Vij is the distance between spins i and j, Zij is the projection of the corresponding 
displacement onto the 2-axis, the easy axis of the spins. Finally, cxj is an Ising spin variable 
such that (7j = ±1 when the true spin on site i is ±5*. 

Since the dipole field is long ranged, and Edm ^ Edn, the fiip of the ith spin changes the 
bias field on a large number of neighbouring spins, and thus changes the fiip probability for 
those spins significantly. The relaxation of the magnetization of the entire solid is therefore a 
complex coupled process in which every individual spin essentially waits until it experiences 
a bias field less than W in magnitude, and then fiips with a probability per unit time 



equal to approximately jW . The flip of this spin changes the bias fleld at many other 
molecules, and if one of them then happens to have a near- zero bias fleld, it flips, leading to 
the possibility of flips at yet more molecules. Ref. j4| refers to this scenario as a long-range 
Glauber model. 



III. MONTE CARLO SIMULATION 



A. Simulation protocol 

As explained in Sec. HI in this paper we only report on simulations on spherical samples of 

spins on a cubic lattice in order to eliminate the effects of inhomogeneous demagnetizing 
flelds. In addition, we only consider the demagnetization process. Thus, the spin Oi is 
initialized to the value +1 at every site. Starting from this conflguration, we simulate the 
time evolution of the sample (as described below) for between 60 and 500 runs, and then 
average the total magnetization of the entire sample over these runs. We have performed 
simulations for two sample sizes, with = 9, 171 and 82,519. 

The initial spin polarization creates an almost delta-function-like distribution of bias 
flelds centered at zero fleld, exactly as expected theoretically. We see small deviations from 
a perfectly uniform distribution due to the flnite size of the sphere. 

The evolution of the system from time step t to the next time step t + is carried out 
using the following protocol. At time t, the bias energy Ei is computed at every site using 
Eq. (12. 3p . All spins are then flipped or not flipped using the flipping protocol described 
below. We are now at time t + dt. The bias flelds are recomputed at all sites, and the 
process is repeated. 

The flipping protocol we employ entails a slightly modifled flip probability 

Pmp. = ^e(H^-|E,|)dt, (3.1) 

instead of the original form (12. 2p . Here, 0(-) is the Heavyside function equal to unity for 
positive argument and to zero for negative argument. In other words, a spin flips only if 
the bias fleld on it is less than W in magnitude. This modiflcation is not material to the 
physics, and it reduces the run time of the simulations. We refer to the spins in the window 
\Ei\ < W as reversible. We have also used Eq. (12.21) in a few cases, and not found any 
signiflcant differences in the results. Further, A2 = i/vr^' ^^'^ prefactor in Eq. (13. ip 
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is chosen to ensure that the integral J^^piiip{E)dE is unchanged. In this way, the total 
magnetization that flips in a large subvolume containing many spins is unaffected. For 
future use we deflne 

ro = ||- (3-2) 

An important consideration arises with regard to the values of Eam, A, and W to be used 
in the simulation. We know that the ratio of these quantities for real Feg is Edm/ ^ ~ 10^, 
and Edm/W ~ 10. Due to the long ranged nature of the dipole fleld, when a spin flips, it 
has the potential to bring ~ lOEdm/W spins into the reversibility region \Ei\ < W. We 
refer to this as the influence sphere of the spin. To overcome finite size effects, we must 
make sure that our simulation includes a large number of inffuence spheres. Secondly, the 
rate at which a spin ffips, even if it is within the reversibility window, is governed by A, 
and our simulation would be much too slow if we used the actual value of A/ Edm- We 
have therefore chosen different values for these quantities while still ensuring the physically 
important restriction Edm ^ W ^ A. Specifically, we take A2 = 2.0, Edm = 5OA2, and 
vary W over a range of values between A2 and Edm- 

The next consideration is over the choice of the time step dt. We set dt — O.OlEdm/ 
and hence independent of W. This is done in order to remain true to the idea that the flip 
probability for a reversible spin should depend on W only through the rate Fq, and not dt. 
With our choice of dt this probabihty is 

Pfiip = ^dt 

= 01 ^2 Edm 

Uw Ai 

= 0.01:^. (3.3) 

By choosing Edm/4W < 10, we ensure that the flip probabihty in one time step is not 
too large, which in turn ensures that our discretization of time is not too coarse, and that 
the simulation is sufficiently close to a continuous process. At the same time, pfnp is large 
enough that we do not expend unnecessary time steps in waiting for the spin configuration 
to change by a meaningful amount. The time-scale r = Edm/ ^2 demarcates short versus 
long times, and we shall study relaxation for ~ lO^r in some cases, i.e., ~ 10^ time steps. 
For real Fes we have r ~ 10^ sees. 

Some other details of the simulation are as follows. The spherical sample is built from a 
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cube having an odd number of sites on a simple cubic lattice with lattice constant 'a', and 
selecting those sites within a distance Da/ 2 of the origin in order to get a sphere of diameter 
Da. The two system sizes = 9,171 and = 82,519 correspond to sphere diameters 
D = 27 and D = 55 respectively. The sites are indexed from 1 to A^, and their Cartesian 
coordinates are stored in one- dimensional arrays. To reduce computer time, at the start a 
one- dimensional look-up table is made of the kernel Kij by converting the triple of distances 
{xij, Hij, Zij) into a single unique number using some artificial but easy-to- implement formula 
that is invertible, i.e., capable of yielding the triple {xij,yij, Zij) from the single number. 

B. Quantities Measured 

The central quantity of interest that is measured in our simulations is the magnetization, 

m={N^- N^)/N, (3.4) 

where A'^ and A^^ are the number of up and down spins. The magnetization is measured at 
every time step. 

In addition, we also measure at every time step, the bias distribution p{E), defined 
such that p{E)dE is the fraction of spins experiencing a bias field between E and E + dE. 
The bin width for numerical purposes is chosen as W itself as this is a sufficiently small 
number compared to Edm- Secondly, the distribution is measured for biases that satisfy 
\Ei\ < 15Edm- In practice, we find that the fraction of sites that lie outside this range is 
0(10-2). 

C. Results of the Simulations 

As mentioned above, we have performed the simulations for different relative values of 
W and Edm- For a test case, we made the contraphysical choice W ^ Edm- In this case 
we expect each spin to remain reversible most of the time, and rarely move out of the 
reversibility window when neighbouring spins flip. Each spin should then relax essentially 
independently of the others, leading to exponential relaxation of the magnetization with a 
rate 2ro. This is indeed what is observed, giving us confidence in our numerical code. 

The physically interesting simulations are performed for Edm 3> W. In Fig. [H we show 
the magnetization versus time for one such simulation over a time lOOOr. It is evident that 
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the decay of m is nonexponential, and that there is a steep initial drop in m over a time of 
order r. This drop is shown in more detail in Fig. [21 and is quite well fit by a square-root 
form; we discuss this in more detail in Sec. I VII In both these figures, we have performed an 
average over 60 runs. 

In Figs. 13] and S] we show the short- and long-time bias distribution p{E) for the same 
parameters as in Figs. [1] and [2l At short times, the distribution is marked by three clear 
peaks, as well as a few shoulders, which we shall explain in more detail in Sec. |Vl Here 
we note that the two main peaks other than at the center are at —AEdm and 8Edm- It 
is also evident that the peaks and shoulders become less distinct as t increases. Indeed, 
for t > lOOr, they disappear completely, as shown in Fig. SJ Here we see a new feature 
developing, namely a hole in the distribution at = 0, for t > 500r. 

The bias distributions also provide a good indicator of whether our system size is large 
enough and whether the averaging procedure is valid. To this end, we show in Figs. [5] and 
|6] the short- and long-time distributions for the smaller sample size (A^ = 9171) but all 
parameters the same as in previous figures. The two figures are drawn for averages over 60 
and 30 runs, respectively. As can be seen, the statistical scatter is only minimally greater, 
and the quantitative features — heights and locations of the peaks at short times, the hole at 
zero bias at long times — are identical. Finally, in Fig. [TJ we show the short-time distribution 
for a single run of the larger sample. The features seen in the 60-run average are all clearly 
present, showing that questions of self-averaging do not arise in this system. 

IV. RATE EQUATIONS FOR MAGNETIZATION RELAXATION 

To understand our simulations, we have developed a theory based on rate equations. The 
key realization lies in the very different role played by the reversible and the nonreversible 
spins, and that we therefore need to understand the time-development of each set separately. 
We denote by N^, A^t; and A^j, the total number of reversible spins at any instant (i.e., those 
with a bias satisfying l-E] < W), and the parts of this number whose spins are up or down. 
Corresponding lower case symbols n^, n^^, etc. are used for the fractions Nr/N, Nr^/N, 
N^/N, etc. We also denote the number of nonreversible spins, A^ — A^^, by Nf, and the 
sets of spins of various types by Sr, Sr^, Sp etc. These sets obey obvious relations such as 
Sr = 5r.-|- U Sri ^^'^ ^'^ which need not be listed. It also pays to introduce the reversible 
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magnetization, 

rur = n^t - Url, 

the total magnetic moment M = Nm, and its reversible part, = Nrrir. 



(4.1) 



A. Processes that change the state of a spin 

We now examine how different spins can develop in a small time interval dt. A non- 
reversible spin (at site i, say) can 

1. Move into the reversible bias range with a probability pin,i. 

2. Remain in the non- reversible range with a probability 1 — Pin,i- 

Naturally, since this spin cannot flip in the interval dt, these possibilities depend on 
the behavior of other spins. We shall address the probability pin,i below. 

A reversible spin (again taken to be at site i), on the other hand, can do the following: 

1. Flip and move out of the reversible range with probability PflipPout,i- 

2. Flip and remain in the reversible range with probability Pmp{^ — Povit,i)- 

3. Not flip and become nonreversible with probability (1 — Pflip)Pout,i- 

4. Not flip and stay reversible with reversibility (1 — Pflip)(l — Pout.i)- 

Once again, the probability Pont.i depends on the behavior of other spins, and will be 
estimated below. We have also introduced the quantity 

Pflip = Todt, (4.2) 

in which the index i is omitted in pf^p, since this is the flip probability for all reversible 
spins. Clearly, our model assumes that the processes of flipping and of moving in or out of 
the reversibility range are independent, which in turn means that different spins flip or do 
not flip completely independently of each other, with a probability that depends only on the 
local bias. This assumption will be valid provided the bias distribution p{E) is reasonably 
spatially homegeneous across the sample at all times. Such is the case for our spherical 
samples, but will need to be reexamined for other shapes. 
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With the above hypothesis, the change in the numbers of various types of spins in a short 
time interval dt are easily written down. For dNrf, we have 

dNrf = - [pflipPout,i+Pfiip(l-Pout,i) + (l-Pfiip)Pout,i]+ Pfiip(l-Pout,i)+ Y Pi^,i- (^-S) 

The first four terms on the right correspond to the four processes enumerated above for re- 
versible spins, while the fifth term corresponds to nonreversible up spins becoming reversible. 
Simplifying, we get 

dNr^ = - Y [PfliP + (1 - Pflip)Pout,j + Y Pflip(l -Pout.i) + Y Pin.i- (4-4) 

Similarly, 

dNrl = - Y [PfliP + (1 -PRip)Pout,i + Y Pflip(l " Pout.i) + Y Pin.i- (4-5) 

Adding the last two equations, we get a very simple equation for the change in the total 
number of reversible spins, 

dNr = -Y Pont,i + Pin^i^ (4-6) 

which does not depend on pflip at all, since we do not discriminate between up and down 
spins in the set Sr, and the changes in its size are a function of the behavior of neighboring 
spins of the members of this set. 

By taking the difference of Eqs. (14.41) and (14.51) . we get the change in the unnormalized 
reversible magnetization: 

dMr = dNr^ - dNrl- (4.7) 

We can simplify the expression that results upon substitution of the actual forms of dNr^^ 
and dNrl anticipating that the probabilities Pout,i and pin.i will also be proportional to 
dt. Thus terms such as PflipPout.i are 0{dtY and may be omitted. In this way, we get 

dMr = -2pflipM,, - {Y Po^i^i ~ Pout,i) + {Y Pin.i " Y P^^^i) ' (4-8) 

iGSr^ iGSr^ i£Sr^ IGS^^ 

Lastly, we find dM, the change in the total unnormalized magnetization. Since this 
change can come about only by the fiipping of reversible spins, and since each fiip changes 
M by 2, 

dM = -2Y Pflip + 2 Y PfliP 

= -2pflipM,. (4.9) 
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B. The probabilites and pout 



For the equations for dNr, dMr, and dM to be useful, we need the probabihties pm and 
Pout- Let us begin by considering a nonreversible spin i G Sr that sees a bias Ei > W. For 
this spin to move into the reversible range, reversible spins at other sites will need to flip 
and alter the bias at site i to satisfy \E[\ < W, where the prime indicates the bias after a 
time interval dt. Now, 

E[ = Ei + J2' K^JdaJ. (4.10) 

Here, the site i is excluded from the sum, and daj is the change in the spin at site j in 
the time dt. The requirement that \E[\ < W implies that only a particular set of reversible 
spins determined by the geometry of the lattice and the form of the dipole kernel Kij can 
be effective in making spin i reversible. We shall refer to such spins as triggering spins. To 
estimate their number we make the critical simplification that we may ignore simultaneous 
spin flips since such processes will have a very low probability proportional to (c/t)^, which 
may be neglected as dt is infinitesimal. Thus, in Eq. f l4.10p . we take dcfj = for all but one 
distant reversible spin. Taking this spin to be up, so that daj = —2, we get 

E'^ = E,-2K,,i^), (4.11) 

where the arrow in -ft'jj(t) indicates that the distant spin is up. The condition \E^\ < W 
then implies that 

^ < < (4.12, 

Similarly, if the distant spin is down, we require 

< K,M S (4.13) 

We now find the number of sites for which the couplings Kij lie in the range fl4.12p or 
dHS]). If we define 

K,{E,) = ^{E,-W), K2iE,) = ^{E, + W), (4.14) 

then these two ranges correspond to intervals [Ki,K2], and [—K2, ~Ki\ in which Kij must 
lie. Let us denote the numbers of sites in each interval by N^Ki,K2] ^-nd N[_k2-Ki]- We have 

N[K„K2] = / 9iK)dK, (4.15) 
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and similarly for Ny_K2-Ki\i where g{K) is the density of dipole couplings found in Ref. |13 |. 
That is, g{K)dK is the number of sites for which the coupling to a central site lies between 
K and K + dK. We have 

g(K)-a^, a = ^. (4.16) 

It then follows that 

= "^""ii - it) 

Since the number of distant sites at which a triggering spin could be located is independent 
of whether that spin is up or down, we can calculate the probability that spin i will become 
reversible, that is to say Pin,i, as the product of three factors: (i) the number (14.181) . (ii) the 
fraction of these sites at which the spin is itself reversible, n^, and (iii) the probability that 
any one of these spins will flip, pfup. Thus, 

Pin^i = AaUr^^^^Todt. (4.19) 

The above calculation assumes once again that the local reversible fraction rij. in the vicinity 
of spin i is spatially homogeneous, and thus independent of the location of site i. 

We next turn to the calculation of Pout, which proceeds in close parallel to that of pi^- 
Consider a reversible spin at site i, i.e., the bias Ei obeys \Ei\ < W. We refer to this site 
as the central reversible spin. This spin will become nonreversible if a distant reversible 
spin flips in such a way as to push the bias at site i outside the interval [—W, W]. Suppose 
the distant spin flips from up to down. Since we have \Ei\ < W and want \E'j^\ > W, the 
coupling Kij must be such that 

K,,^[K,,K2], (4.20) 

where Ki and K2 defined in Eq. (l4rT4|) . Noting that now Ki < and K2 > 0, the 

number of sites that meet this requirement is given by 

. E,.mW 

IK2 

Similarly, if the distant spin flips from down to up, the condition on Kij is 



K,j^[-K2,-K,], (4.22) 
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which is met by a number of sites equal to 

^ g{K) dK + 1^ g{K) dK = 4a:^^^|^, (4.23) 



which is the same as Eq. (14.211) . Thus, once again, the number of sites on which a triggering 
spin can be located is independent of whether that spin is up or down, and we may calculate 
Pout,?, as the product of (i) the number of sites (I4.2ip . (ii) the fraction rir that the spin on 
one of these sites is reversible, and (iii) the probabibility Todt that this spin will indeed flip. 
Thus, 

Pout,i = 4Q;nr. ^'^'^^^ ro(jt. (4.24) 

The expressions (I4.19P and (I4.24p suffer from unpleasant singularities when Ei = ±W. 
These singularities are unphysical, and are a consequence of using the modified spin-flip 
probabillity (13. ip with the hard cutoffs at ±W. Better estimates are obtained by noting 
that for Pin,i, \Ei\ is likely to be much bigger than W, while for Pout,i the converse is true. We 
therefore neglect the term in the denominator of Eq. fl4.19p and Ef in the denominator 
of Eq. (I4.24p . leading to the expressions 

pUEi) = Aarir^^Todt, (4.25) 

i 

PoutiEi) = AaUr^Todt. (4.26) 

We note here the intuitively reasonable fact that pout is much greater than pi^. The set Sr is 
much smaller than Sf, so an initially reversible spin will be knocked out of reversibility by 
almost all flips of neighboring spins. By contrast, to move an initially nonreversible spin into 
reversibility, one must cancel the preexisting bias at the nonreversible site nearly exactly, 
which can only be done by flipping distant spins at a very specific set of sites. For this same 
reason, Pout,i essentially does not depend on E^, while pin^j does. 

Note also that and pout are both proportional to dt as anticipated earlier. 



C. The rate equations 

We now substitute Eqs. (1419|) and (^M> into Eqs. (gSD, dM]), and (gl]) for dN^, dM^, 
and dM, and divide by the total number of spins at the same time in order to get equations 
for intensive quantities. Let us begin by considering the two sums in Eq. (14. 6 p one by one. 
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Since Pout,i is independent of Ei as noted above, we have 

^ Pont,i = Pout^ = 4:anl^^Todt. (4.27) 

For the second sum, we need to sum over the set Sf- We do this by including all sites where 
the bias exceeds W in magnitude. This leads to the approximation 

1 E Pin,. = 'iaUr^J^Todt, (4.28) 

where J-" is a dimensionless functional of the bias distribution p{E), given by 

J^[p{E)] =W' [ P^dE. (4.29) 

J\E\>W E^ 

Hence, 

dUr Edm , ^ .X 

— ^ = -4«ro-^n,K - ^ • 4.30 
dt W 

Next, we examine Eq. (14. 8 p for dMr/dt. For the term with the sums over the sets Sr^ 
and Sr^, we have, 

^( H Pont,i - J2 Pout,i) = 4:anrmr^^Todt. (4.31) 

For the remaining two sums, we estimate the sizes of the sets Sp^ and Spi as and 
times the size of Sf on the theory that when <^ 1, most of the spins are nonreversible 
and the bias at any site is uncorrelated with whether the spin at that site is up or down, 
and that when — 1, m ~ m^. It follows that 

1 / ^ - \ . Edm 



^ ~ ^ ^''^'V ~ 4:anrm-^J'Todt. (4.32) 
Hence, 

= —2TQmr — AaTQ-^nr{mr — nij^). (4.33) 
Lastly, we obtain the equation for dm/dt, which is the simplest of all: 

dm ^ , . 

— = -2Tomr. 4.34 
dt 

Equations fl4.30p . f l4.33p . and f l4.34p are the desired rate equations. They are manifestly 
nonlinear, but more importantly and contrary to our initial hope, they are not a closed 
system because of the presence of the functional J-" of the full bias distribution p{E). At 
present this puts a big limitation on their use. For the relaxation problem we have been 
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able to circumvent this limitation by constructing an interpolation form for p{E) which we 
believe is reasonably accurate and self-consistent over a wide range of times, well past that 
over which the square-root time development is seen. We describe our approximation for 
p{E) in the next section. 

V. THE BIAS DISTRIBUTION 

A. The three- Gaussian approximation 

As seen from the Monte Carlo simulations, the bias distribution at short times is domi- 
nated by three peaks at £■ = 0, £■ = —AEdm, and E — SEdm- The locations of the two side 
peaks are a strong indicator of their origin. Consider a site with its six nearest neighbors. 
Four of these neighbours are in the xy plane, and two are along the z axis. If any of the 
neighboring spins in the xy plane flips from up to down, the bias at the central site will 
change by an amount —AEdm, while if any of the z axis neighbours flips, the field at the 
central site will change by SEdm- This explains the peak locations. Further, since there are 
twice as many near neighbours of any site in the xy plane as there are along the z axis, we 
should expext the peak at —AEdm to be about twice as high as the peak at SE'^^ as long as 
Nr <^ N. This is also seen in the data. The smaller peak at —8Eam and shoulder at 4£'ci^ 
can also be associated with spin fiips at pairs of near neighbour sites. 

Motivated by this idea, we try and represent p{E) as a sum of three Gaussians centered 
at 0, —4Edm, and SEdm- Suppose that at a given time, N_i_ spins have flipped where A^^ -C N, 
allowing us to ignore the possibility that two flipped spins are near neighbours of each other 
or even of a common third spin. Then there are AN^ spins that have a flipped neignbour in 
the xy plane, and 2Ni spins that have a flipped neigbour along the z axis, leaving N — GN^ 
spins which have no fiipped neighbours at all. Thus the weights of the 0, —AEdm and SEdm 
peaks are proportional to (1 — Gn^,), 471^,, and 2ni respectively. We can further argue that 

1/2 

the widths of all three peaks are equal and proportional to , since the fields at sites far 
away from all flipped spins should continue to vanish on average, but should have a variance 
that grows linearly with the number of flipped spins. For a site next to a flipped spin, this 
variance is simply realized around the shift produced by the flipped neighbour. Thus for 
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rii <^ 1, the three- Gaussian approximation (TGA) to p{E) takes the form 

p{E) ^ (1 - 6n^)go{E) + 2n^g+{E) + An^g_{E), (5.1) 

where (with a = 0, +, or — , and Eq = 0, E^ = SE^m, and E_ = —AEdm) 

^„(E) = (27^n,a2)-l/2e-(^-^")^/2'^^^^ (n^ « 1). (5.2) 

The quantity a is Edm times an unknown constant of order unity. 

The arguments underlying Eq. (15. ip start to become questionable for as small as 0.1, 
since sites with two near neighbour flipped spins start to become significant. To enable us 
to consider larger values of rii, we generalize the TGA to the form 

p{E) ~ aogo{E) + a+g+{E) + a.g.{E), (5.3) 

where 

go^{E) = (27ra2)-i/2e-(^-£^^)V2.^_ (54) 

That is, the peaks of the three Gaussians are still taken to be at 0, —AE^n and SE^n, the 

1/2 

widths are taken to have a common value a not necessarily proportional to , and the 
weights Oo, a+, and a_ are allowed to become arbitrary. We will determine these weights 
and the width by the procedure described in the next subsection. The form (15.11) at small 
ni will serve as a check on the procedure. 

It is apparent that the TGA is qualitatively incapable of accounting for the very narrow 
hole that is burned in the distribution at long times, but here a different approximation 
scheme can be developed as the origin of the hole is physically obvious. 

B. Moments of the bias distribution for uncorelated spins 

Our discussion above implies that for very small ni, the flipped spins are randomly 
distributed in the lattice without any spatial correlations. We therefore extend this idea to 
larger ni and consider a model in which the spin on each site is up or down independently of 
other spins, with probabilities (1 ± m)/2, where m is the magnetization. We then calculate 
the first three moments of this model, and match those to the moments of the TGA, Eq. (15.31) . 
These three moments, plus the normalization (or zeroth moment) give us the four conditions 
needed to determine the four quantities and a. 

16 



The bias at any site i is given by 

E, = J2K,,a,. (5.5) 

Consider first the uniform spin configuration with m = 1, i.e., cTj = 1 for all i. We know 
that in this case the bias vanishes at all sites except those in a narrow layer near the surface 
of our spherical sample. Hence we may take 

E^u = (5-6) 

for essentially all sites. This result will be employed repeatedly in the calculations of the 
moments for configurations in which m ^ 1. Thus, for the first moment, we have 

= 0. (5.7) 
Similarly, for the second moment, we get 

{Ef) = J2' K,,K,k{cJjCrk). (5.8) 

The prime on the sum signifies that j ^ i and k ^ i. Now equals 1 if j = k, and 

if j ^ k. Hence, 

= '^2EUl-m'), (5.9) 
where we have used Eq. (15. 6p . and defined 

'^2 = T^E^5- (5-10) 
Numerical evaluation of the sum gives 

K2 = 53.427 (5.11) 
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For the third moment, we have 



{Ef) = Y! KijKikKu{(Tj(rkai). 



(5.12) 



Again, the prime signifies that j ^ i, k ^ i, and / 7^ i. The only issue requiring care in 
performing the sum is the enumeration of the various cases of equahty or inequahty of the 
indices j, k, and I. The first case is where all three indices are distinct. Then (ajakcri) = m^, 
and the contribution of this case to {Ef) can be evaluated as 



m 



m 



m 



(E' ' - 3 E' E' + 3 E' - E' ^' 



2m^ J2 K. 



(5.13) 



In line 2 above we have used the symmetry of the summand, and in line 4 we have used 
Eq. dSH). 

The second case is where two of the indices j, k, and I are the same, but distinct from 
the third. Now {(Jjcrkcri) = m. This case has three identically contributing subcases, and for 
its net contribution to (Ef) we have 



'ji) 



3m 



E'^E'^^-E'^ 



-3m J2 



(5.14) 



where we have again used Eq. (15. 6p in the last line. 

The third and last case is that where j = k = I. Now {ajatai) = m, and the contribution 
to {Ef) is, therefore, 

{Ef), = mJ2Kf^. (5.15) 
Adding together Eqs. f l5.13p . f l5.14p . and f l5.15p . we get 



{Ef) = -2m{l-m')J2Kf^. 



(5.16) 
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We write this as 

{Ef)=KsEl^m{l-m'), (5.17) 

where 

«^3 = -J-E^5 = 190.47, (5.18) 

and the last resuh is found numerically. 

It should be noted that in this model, the moments of E are simply geometrical constants 
determined by the type of lattice times the appropriate power of the energy scale Edm- 



C. Moment matching 

We now match the moments from the previous subsection with those of the three- Gaussian 
approximation (15. 3p . The latter yields 

(E) = 8Edn.a+-4Edma-, (5.19) 

{E') = a'iao + a+ + a_) + 64^^^+ + IQEj^a., (5.20) 

{E') = l2Ed^a\2a+ - a_) + 512^^^+ - 64ELa- (5.21) 

Equating these moments to those from the uncorrelated spin distribution yields 

4Ed„(2a+ - a_) = 0, (5.22) 

a\a^ + a+ + a_) + 16EL(4a+ + a_) = K^E^^il - m^), (5.23) 

12Eama\2a+ - a_) + 64E3^(8a+ - a_) = K3Ej^m(l - m^). (5.24) 

Solving these equations along with the normalization condition, 

ao + a+ + a_ = 1, (5.25) 



we obtain 



ao = 1- Y^m(l-m2), (5.26) 

a+ = ^m(l-m2), (5.27) 

a- = ^m(l-m2), (5.28) 

a' = ^(4/^2 - ^3m){l - m')Ej^. (5.29) 
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At this point let us ask whether the solution fl5.26p - fl5.29l) approaches Eqs. (15. ip and (15. 2p 
when rij, ^ 1. In that limit, since m = 1 — 2n^, m(l — m^) ^ (1 — m^) = 4n^. Feeding 
in the value K3 = 190.2, we get a+ = 1.98n|, and a„ = 3.96n|, instead of 2n^ and 471^^. 
The differences are rather small, however, and can be eliminated entirely if we make the 
replacement 

Kg ^ 4 = 192. (5.30) 

This leads to the final forms we shall use in our three- Gaussian approximation, Eqs. (15. 3p 
and dEU: 

ao = 1 - |m(l -m^), (5.31) 
a+ = im(l-m2), (5.32) 
a_ = m(l — m^), (5.33) 
a' = (K2-48m)(l-m2)EL- (5-34) 

D. Comparison with simulations 

When we now compare the TGA with the simulations, we discover that the agreement is 
off by ~ 10% if we use the value K2 = 53.4. This value was calculated for an infinite lattice, 
and for a finite sized sample the variance of Ef should be smaller. Using the value 50 
appropriate to the 82519 spin sample, we find that the agreement is considerably improved. 
In Fig. [8] we show the TGA with the choice k,2 = 50 along with the results of the simulations 
for the 82519 spin sample for t/r = 0.1, 0.3, and 0.5, where r = Edm/^"^- At these three 
times, m = 0.93, 0.89, and 0.86. The agreement becomes poorer for larger t, and it is about 
as good as could be expected given how simple-minded the approximation is. 

VI. SHORT-TIME DECAY OF MAGNETIZATION: THE Vi LAW 

In Fig. [2] we show m{t) for short times from our simulations, and from solving the rate 
equations with the value K2 = 53.4. As can be seen the general trend is the same, although 
the detailed agreement is only good to about 3%. Once again, the agreement is improved if 
we set K2 = 50, as shown in Fig. |9l The same data are shown on a log-log plot in Fig. [TOl 
As can be seen, both the simulations and the rate equation show a power law behavior, with 
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the same exponent. The best fit gives an exponent of 0.46, which is very close to 0.5 as it 
would be for -y/t behavior. We now show that this behavior can be understood analytically 
on the basis of our rate equations, and that this exponent does not depend on the choice of 

The first key point is that starting from a delta-function at t = 0, the bias distribution 
becomes broader than the reversibilty region at some ultra-short time when the fraction of 
flipped spins is still very small. From Eq. fl5.29p . we find that for ^ 1, 

^ A'Ej^n^, (6.1) 

where = 4^2 — 1^3 ■ Thus a < W only as long as < {W/AEdmY, which is of order 
10"'^. For such ultrasmall values of n^, = 1 — Qn^, and J-" ~ 0, so the rate equation for 
simplifies to 

This has the solution = (2aro-Edm/3W^)t, and so the condition that a < W holds only 



W 



for t < tus, where 

is an ultra-short time scale of order 10"^r. 

It follows that there is a large range of times, tus ^ tr, for which a ^ W even though 
ni <^ 1, i.e., very few spins are flipped. Thus almost all the weight in the bias distribution is 
still in the central Gaussian, i.e., ~= 1, and the dimensionless functional that determines 
the repopulation of the reversibility region can be approximated as 

TT/2 ^00 p-Ey2a^ 

7 = 2^_ / ^——dE. (6.4) 
Now, by integrating by parts, we get 



Jw E^ W Jw 



J_g-W2/2(72 _ J_ 



The last expression can be expanded in powers of W, and we get 



(6.5) 



2W / [WW , , , 
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The second key point is that even though n^. <C 1, almost all the spins have been knocked 
out of the reversibility region, i.e., ra^ ^ 1. To see this we again approximate the bias 
distribution by neglecting the weight outside the central Gaussian, and setting ao = 1, so 

rir ~ -^L= r e-^'/^'^'dE. (6.7) 

Expanding the integrand in powers of E and integrating, we get 



Thus, to first order in W/a, rir = J-", and the difference is of higher order: 

n,-^=— . (6.9) 
We can express this in terms of rir itself by using Eq. f l6.8p . We have 

— ^ ]l^nr, (6.10) 



so 



rir — T = ^^r- (6-11) 



The rate equation for then reads 



inl, (6.12) 



di 2 

where we have defined 



C = 4aro^ = a^. (6.13) 

The integration of Eq. fl6.12p is elementary. Since this equation only holds for t > tus? we 
can write the integral in the form 

^ = 7rC(t + r), (6.14) 

where t* is a time of order tus- We thus have an explicit solution for the time dependence of 
the reversible fraction: 

The other rate equations can now be solved as follows. We have by definition, 

rir = + Url, TTlr = Tlr'^ — Tlr]^- (6.16) 
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Since Uri < ^ 1, the answer for is immediate: 

1 1 



The rate equation for m now reads 

dm 2To 1 



(6.17) 



(6.18) 



dt v^{t + t*y/^' 

The integration is again elementary. Assuming that tus ^ ^ ^ t, we can write the result as 



m{t) ~ 1 - ^Ti/2t, (6.19) 

where 

= 16^ = (6.20) 

Equation (I6.19P is the experimentally observed \/i law. 

As noted in Sec. [U a very pretty heuristic argument for this result is given in Ref. 
These authors reach the same conclusion by arguing that a(t) must be of the order of the 
typical dipole field when the spins start flipping, and thus proportional to a^/i^{t), where 
i{t) is the typical distance between reversed spins. They then note that a? /l'^{t) oc ni{t), 
so that a{t) oc ni{t). They then estimate mr{t) as W/a{t), from which it follows that 
dm/dt ~ l/n^{t)^ and that n^(t) ~ t^^^. As part of this argument, one has that a{t) ~ t^/"^ 
and that nr{t) ~ t^^/^. We find the same behavior for these quantities, but we arrive 
at it in a different (and more difficult!) way since we did not have enough confidence in 
our understanding of the relation between (j{t) and n^it). Instead, we find the delicate 
noncancellation between rir and J-" in order to first find the differential equation obeyed by 
nr{t), and determine that nr{t) ~ t^^/^, after which the equation for m(t) is elementary. The 
agreement with [4j] gives us encouragement that our procedure is correct, and that we can 
use our more detailed rate equations to analyze other experimental protocols in the future. 
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FIG. 1: (Color online) Long-time decay of magnetization for the N = 82519 spin sample, averaged 
over 60 runs. The parameter values are W = 2.5A2, and Ed-ui = 5OA2. 
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FIG. 2: (Color online) Short-time behavior of the magnetization, with the same parameters as 
Fig. [TJ Also shown is the result from numerical solution of the rate equations. 
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FIG. 3: (Color online) Histogram of the short-time bias distribution for the N = 82519 spin sample, 
averaged over 60 runs, with the same parameters as in Figs. [1] and [2j The bin width in the bias is 
5.0. 
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FIG. 4: (Color online) Same as Fig. [3] but for long times, and averaged over 10 runs only. Note 
the reduced scale on the y axis. 
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FIG. 5: (Color online) Same as Fig. [3] for the = 9171 spin sample. The average is over 60 runs. 
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FIG. 6: (Color online) Same as Fig. |4]for the N = 9171 sample. The average is over 30 runs. 
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FIG. 7: (Color online) Single-run short-time bias distribution for the = 82519 sample, with the 
same parameters as before. There is no averaging. 
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FIG. 8: (Color online) Comparison between the three-Gaussian approximation (TGA) to the bias 
distribution and the simulation results for short times. The sample has N = 82519 spins, and all 
other parameters are as in previous figures. 
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